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ABSTRACT 



We use a block Lanczos algorithm for computing a few of the smallest 
eigenvalues and the corresponding eigenvectors of a large symmetric matrix rather 
than computing all the eigenvalue-eigenvector pairs. The basic Lanczos algorithm 
generates a similar matrix which is block tridiagonal from a given large symmetric 
matrix. The size of the generated tridiagonal matrix depends upon the number of 
the smallest eigenvalues to be computed. The result is savings in computations and 
storage. The block Lanczos algorithm is well-suited for problems involving multiple 
eigenvalues. 

In this thesis, we develop the block Lanczos algorithm to estimate the 
direction-of-arrival (DOA) of a point source based on the observations measured at 
a linear array of sensors and compare the performance with that of a single vector 
Lanczos algorithm. The results of the computer simulation experiments conducted 
with this method are presented and discussed. 
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I. INTRODUCTION 



It is necessary to estimate the power spectral density (PSD) function in 
several signal processing problems. Of particular interest is the estimation of 
spectral lines in noise. One related problem is the direction— of— arrival estimation of 
point sources based on array measurements. Methods based on subspace estimation 
have been proposed in this regard [Ref. 9]. Much of the basic literature on these 
methods has been borrowed from functional and numerical analysis [Refs. 8, 21]. 
Determination of the subspaces needs a complete eigendecomposition of the given 
data autocorrelation matrix. The general eigendecomposition of an n*n matrix 
requires 0(n 3 ) computations. Often we only need to compute a few of either the 
smallest or the largest eigenvalues and the corresponding eigenvectors of a large 
symmetric autocorrelation matrix rather than all the eigenpairs of the matrix. In 
this research, we develop an algorithm for computing a few of the smallest 
eigenvalues and the corresponding eigenvectors and apply this to high resolution 
spectral estimations problems. 

The algorithm we develop is an extension of the method of minimized 
iterations due to Lanczos [Ref. 9]. Paige [Ref. 4] experimented with the Lanczos 
algorithm and found that a few of the extreme eigenvalues of a tridiagonal matrix 
would often converge rapidly to the similar eigenvalues of a real symmetric matrix 
R much before the entire set of eigenvalues are computed. 

The block Lanczos algorithm is an extension of the basic single vector Lanczos 
algorithm. We iterate with a block of vectors rather than with a single vector, and 
generate a reduced block tridiagonal matrix that is similar to the original 
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autocorrelation matrix R. The basic idea is that the eigenvalues of the block 
tridiagonal matrix are approximately the same as the extreme eigenvalues of R. 
Several researchers have reported the block Lanczos algorithm and its variants, in 
particular, Cullum and Willoughby [Ref. 1], Kahan and Parlett [Ref. 6] and Golub 
and Underwood [Ref. 4]. 

The objective of this thesis is to develop a block Lanczos algorithm to extract 
a few of the smallest eigenvalues and then estimate the corresponding eigenvectors. 
The smallest eigenvalues are said to correspond to the noise subspace of the spectral 
or array measurements. The proposed algorithm will be used to estimate the 
spectral lines which may represent the direction-of-arrival of point sources in low 
signal-to-noise ratio environments. The block Lanczos algorithm will have to be 
compared with the single vector case with respect to the spectral line estimation 
performance. With these goals in mind, we now proceed to discuss how the thesis 
is organized. 

In Chapter II we introduce and summarize the single vector Lanczos 
algorithm. Complete and selective reorthogonalization of Lanczos vectors is 
discussed there. In Chapter III, we describe and develop the block Lanczos 
algorithm and present the results of the experiments carried out with a computer 
program implementing this algorithm. Also, in this chapter, we compare the 
performance of DOA estimation of the block Lanczos algorithm with that of the 
single vector Lanczos algorithm. Finally, in the last chapter, we discuss and 
summarize the results of the Lanczos method and also make some recommendations 
for future work. 
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II. SINGLE VECTOR LANCZOS ALGORITHM 



The single vector Lanczos algorithm is used to tridiagonalize a real symmetric 
matrix R. This algorithm is based on the concepts such as Krylov sequences and 
subspaces, orthogonal projections of matrices, and Ritz vectors [Ref. 1]. In this 
chapter we describe the general Lanczos recursion and a practical Lanczos 
algorithm. Issues related to the reorthogonalization of Lanczos vectors are also 
discussed. 

A. LANCZOS RECURSION 

Given an n*n real symmetrix matrix R and an arbitrary n*l unit norm vector 
kj, we can obtain a sequence of vectors defining the n dimensional Krylov subspace 
as follows [Ref. 1] 



k 2 = Rk, 
k 3 = Rk 2 = R 2 k 




( 2 - 1 ) 



We can now form the Krylov matrix of rank m as follows 



K„ = [k, k 2 k 3 ■ • • kj 
= |k„ Rk„ R 2 k„ • ■ - , R"-*k,] 



(2-2) 
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and the Krylov subspaces J£ m , m=l,2,- • -,n, are then defined as 

^r m (R,kj) = span {k„ Rk„ R 2 k„ • • •, R^kJ. (2-3) 

We now proceed to obtain a tridiagonal matrix which has the same extreme 
eigenvalues as the given symmetric matrix R. If the columns of a matrix Q m are 
orthonormal, then the matrix 



T„ = QlRQ„ (2—4) 

is an m*m tridiagonal matrix which is an orthogonal projection of R onto the space 
spanned by the columns of Q m . Lanczos [Refs. 1, 9, 14] proposed a recursion to 
generate the columns of the matrix Q m = [qj q 2 • • • qj. The vectors q 1} q 2 , • • *, q m 
are referred to as Lanczos vectors. Thus, for a given real symmetric matrix R, the 
Lanczos recursion produces a tridiagonal matrix T m as a projection of R onto the 
corresponding subspace span{Q m }, spanned by the Lanczos vectors generated. 
These subspaces are in fact Krylov subspaces. [Ref. 1] 

The classical Lanczos recursion for R starts with an n*l initial Lanczos vector 
q, which is randomly generated and normalized. Initially we define /?j= 0 and q 0 E 0. 
We then compute cq, /? i+1 and the Lanczos vectors qj for i= 1 ,2, • • - ,m as follows: 





(2-5) 


/ViQi+i = R di - - Adi-i 


(2-6) 
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(2-7) 



/J].i = ql.iRqi- 

The Lanczos matrix T m is generated by appropriately arranging £*; and /? i+1 , where 
aj are the diagonal elements and /? i+1 are the subdiagonal elements of the tridiagonal 
matrix, so that 



T„(i,i) = o, for l<i<m. 



(2-8) 



and 



T„(i,i+1) = T„(i+l,i) = for l<i<m-l, 



(2-9) 



resulting in 



T = 

m 



Of! P 2 
P 2 @3 
@3 Q '3 



• Pm 

An 



( 2 - 10 ) 



The vectors Ojqj and in Eqn(2-6) are the orthogonal projections of vectors 

Rqj onto the two most recently generated Lanczos vectors qj and q^. Thus, 
updated Lanczos vector q i+1 is computed by orthogonalizing the vector Rqj with 
respect to previously computed Lanczos vectors qj and q 5 _ j. The classical Lanczos 
recursion can be condensed into matrix notation by a group of Lanczos vectors Q m 
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and the Lanczos matrix T m . We obtain the following matrix equation 

= V„ + (2-11) 

where e 5 is the unit vector whose 2 th element is 1 and whose other elements are 0. 
Note that, in this recursion, R is never modified. Also, storage is needed only for 
the Lanczos vectors q iM , q 5 and q i+1 , the Lanczos matrix T m , and the given matrix 
R. 

The classical Lanczos procedure attempts to maintain the orthogonality of the 
Lanczos vectors. Due to the roundoff and other numerical errors, however, 
orthogonality between the Lanczos vectors can only be maintained by incorporating 
some kind of explicit reorthogonalization as the Lanczos vectors are computed. 
Note that the reorthogonalization of these vectors requires extra storage for keeping 
all of the Lanczos vectors as well as additional computations. We discuss the 
reorthogonalization in a later section. 

B. PRACTICAL LANCZOS ALGORITHM 

In this section we focus on a Lanczos algorithm with no explicit 
reorthogonalization incorporated. The loss of orthogonality, if it occurs, may result 
in spurious eigenvalues. Although the orthogonality of the Lanczos vectors is lost, 
some of the eigenvalues of R will still appear as eigenvalues of the Lanczos matrix if 
we make the tridiagonal matrix large enough. Since we are generally interested in 
only a few of the smallest eigenvalues and their corresponding eigenvectors, the loss 
of orthogonality does not critically affect finding approximate eigenvalues and 
corresponding eigenvectors of R [Ref. 1]. Nevertheless, it is possible to find accurate 
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eigenvalues and eigenvectors by using the method in an iterative scheme using no 
reorthogonalization at all, even in the face of total loss of orthogonality [Ref. 12]. 
Now, we present the practical single vector Lanczos algorithm. 

To generate the tridiagonal Lanczos matrix T m we should compute cq and /? i+1 
which are the diagonal and sub-diagonal elements of the tridiagonal matrix, where 
i=l,2, * • - ,m. Given an n*n real symmetric matrix R, we start with an n*l initial 
arbitrary vector qj which is normalized such that I Qi 1 2 = ^ 35 ' n Eqn(2— 1). We now 
define an intermediate vector 



Uj — Rq, (2—12) 

and initialize a o =0 and /? 0 =0. The single vector recursion is then carried out for 
i=l,2,- • • ,m and can be summarized as follows: 



j? 

II 

c 


(2-13) 


w i = Uj - a jqj 


(2-14) 


r i + i = + J 


(2-15) 


Qi + i = w i/ r i + i 


(2-16) 


fo = «iT R< ii + i 


(2-17) 


u i+ i = Rq i+1 — P[<ii 


(2-18) 
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where u ; and Wj are the orthonormal projection vectors for the Lanczos vector qj. 
Here Eqn(2— 13) and Eqn(2— 17) use the modified Gram-Schmidt orthogonalization 
to compute the coefficients a- } and /? i+1 , respectively. If q; is orthogonal to q^, 
then q i+1 will theoretically be orthogonal to both q^ and q 5 . This algorithm is 
quite popular despite the requirement for small amounts of extra computations 
[Ref. 12]. Now the tridiagonal matrix T m is generated by simply filling it with 
and for its entries as shown in Eqn(2— 10). 

1 . Eigenvalue Computation 

To find the eigenvalues ^ of the mxm Lanczos matrix T m , we may use 
the bisection method and Sturm sequencing [Ref. 5]. Actually, one could obtain 
both eigenvalues and eigenvectors of T m by using such methods as the QR algorithm 
[Ref. 16]. However, since we need only a few of the smallest eigenvalues of T m , we 
choose the bisection method to compute them as detailed in the following. 

Given the tridiagonal matrix T m , we define the characteristic 
polynomials p 0 (/z), Pj(//,), • • •, p m (/r) as 



for j=l,2, •••,m. Expanding the determinant in Eqn(2— 20) yields the recursive 
expression [Ref. 5: pp. 305—307]: 



PoM = 1 



(2-19) 



Pj(/r) = det(T m -//I) 



( 2 - 20 ) 



Pj(/0 = (<>j -/^)Pj-i(^) -/7jPj- 2 (^)> j=2,3,- • - ,m. (2-21) 
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The zeros of the polynomial p m (/i) are the eigenvalues of the tridiagonal matrix T m . 
Here we are interested in only a few of the smallest eigenvalues of R. In order to 
compute these values we first define a range [x, y] in which all the desired 
eigenvalues lie. We then carry' out the following iteration to implement the Sturm 
sequencing property. If p(x)p(y) < 0 and x < y, then the iteration 

I x — y | > c( | x | + | y | ) 



y = H if p m (x)p m (y) < 0 (2-22) 

x = P if P m (x)p m (y) > 0 

is guaranteed to converge to a zero of p m (/*), to an eigenvalue of T m . The value 
c is the machine unit roundoff error and the limits on the range [x, y] are given by 



x = min ttj - |/?j| - | /^.jI 
i 

y = max a- + |0j| + I0j.il 



(2-23) 



where we have 0 O = 0 m+1 = 0. The bisection method computes the eigenvalues with 
small relative error, regardless of their magnitude. This is in contrast to the 
tridiagonal QR iteration, where the computed eigenvalues can have only small 
absolute error. 
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2. Eigenvector Computation 

Now, we need to compute the eigenvectors of R to complete the 
eigenpair estimation problem. There are two techniques to compute the 
approximate eigenvectors of R. The first technique uses the Ritz vector Xj and the 
corresponding algorithm uses the relationship [Ref. 14] 

Xj = Q„Zj, j=l,2,---,s (2-24) 

where Zj is one of the eigenvectors of T m , Q m is the matrix of Lanczos vectors, and 

s < -t?-. 

The second technique to obtain the eigenvectors of R corresponding to 
the selected eigenvalues of T m uses the Rayleigh quotient iteration. The iteration 
makes the Rayleigh quotient r(x k ) converge to /q which is one of the smallest 
eigenvalues of T m . The Rayleigh quotient of an eigenvector x k which minimizes 
I (R — r k I)x k || 2 is given by 



T 

x,.Rx k 

r k = r (* k ) = f ' ( 2 “ 25 ) 

x k x k 

As r(x k ) approaches an eigenvalue /q of T m , then the solution to (R — r k I)x k = b k 
will be an approximate eigenvector by using the inverse iteration theory where b k is 
a vector close to zero [Ref. 14]. 

We now briefly present the inverse iteration algorithm. First, we pick 
an arbitrary unit vector x^ then the iteration proceeds as follows 
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